

clear all
*Set version of stata
version 14.2

*seed randomly generated by: https://www.random.org/ (b/t 1 and 1000000000)
set seed 24205253

*Generate datafile to save results:
*PATH NEEDS TO BE UPDATED
cd "..."
tempname pfile
tempfile lrmfile
postfile `pfile' str20(variable) lrm se_lrm t using "`lrmfile'"

*generate 5,000 datasets
forvalues i=1/5000 {
	drop _all
*1st 200 obs will be dropped, T=5,000
quietly set obs 5200
quietly	gen t = _n
quietly	tsset t
*All variables are I(1)
*only x and y are related. x2-x9 are unrelated to y
quietly	gen x = rnormal() if t==1
quietly	replace x=l.x + rnormal() if t>1
quietly	gen x2 = rnormal() if t==1
quietly	replace x2=l.x2 + rnormal() if t>1
quietly	gen x3 = rnormal() if t==1
quietly	replace x3=l.x3 + rnormal() if t>1
quietly	gen x4 = rnormal() if t==1
quietly	replace x4=l.x4 + rnormal() if t>1
quietly	gen x5 = rnormal() if t==1
quietly	replace x5=l.x5 + rnormal() if t>1
quietly	gen x6 = rnormal() if t==1
quietly	replace x6=l.x6 + rnormal() if t>1
quietly	gen x7 = rnormal() if t==1
quietly	replace x7=l.x7 + rnormal() if t>1
quietly	gen x8 = rnormal() if t==1
quietly	replace x8=l.x8 + rnormal() if t>1
quietly	gen x9 = rnormal() if t==1
quietly	replace x9=l.x9 + rnormal() if t>1
quietly	gen g = rnormal() if t==1
quietly	replace g = .6*l.g + rnormal() if t>1
quietly	gen y = x + g

******************************
**Generate lagged variable for each predictor w/ long run effect
******************************
foreach x in y x x2 x3 x4 x5 x6 x7 x8 x9 {
quietly gen l_`x' = l.`x'
}
*************

*drop 1st 200 observations to ensure starting value doesn't
*influence results
quietly drop if t<=200

**Must include lagged values before differenced values, so correct elements of the vce are called
quietly reg d.y l_y l_x l_x2 l_x3 l_x4 l_x5 l_x6 l_x7 l_x8 l_x9 d.x d.x2 d.x3 d.x4 d.x5 d.x6 d.x7 d.x8 d.x9 
matrix vcov =e(V)

*gen local variable of value 1
local i = 1
*loop over variables for which we will calculate a long term effect
foreach x in l_x l_x2 l_x3 l_x4 l_x5 l_x6 l_x7 l_x8 l_x9 {
quietly sum `x' if e(sample)
local lrm = _b[`x']/abs(_b[l_y])
local selrm = sqrt(((1/(_b[l_y]^2))*(_se[`x']^2))+((_b[`x']^2)/(_b[l_y]^4)*(_se[l_y]^2))-(2*(_b[`x']/(_b[l_y]^3))*(vcov[(`i'+1),1])))
*add 1 to i each loop, so appropriate element of vcov is identified.
local i = `i' + 1
post `pfile'  ("`x'") (`lrm') (`selrm') (`lrm'/`selrm') 
}
}


postclose `pfile'
use "`lrmfile'", clear


*identify spurious relationsips for x2-x9
sum t if abs(t)>1.96 & variable != "l_x"
*save # spurious LRMs
scalar sig = r(N)

*identify total simulations for this group
sum t if variable != "l_x"
scalar total = r(N)

*identify percent of spurious LRMs
scalar spurious = 100*sig/total
di spurious

*identify mean lrm for unrelated predictors (x2-x9)
sum lrm if variable != "l_x"


***********************************
***********************************
***********************************
***********************************


